Magnetic Resonance (MR) Imaging has been widely used in medical diagnosis because of its non-invasive manner and excellent depiction of soft tissue changes. Recent developments in compressive sensing theory (see, e.g., Candes, et al., Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory, 52, 489-509 (2006) (Candes et al., 2006, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein), Donoho, Compressed sensing, IEEE Transactions on Information Theory, 52, 1289-1306 (2009) (Donoho, 2006, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein) show that it is possible to accurately reconstruct, e.g. Magnetic Resonance (“MR”) Images from, e.g., highly under-sampled K-space data and therefore significantly reduce the scanning duration.
Suppose x is an MR image and R is a partial Fourier transform, the sampling measurement b of x in K-space can be defined as b=Rx. The compressed MR image reconstruction problem (“problem (1)”) can be viewed as one to reconstruct x given the measurement b and the sampling matrix R. Motivated by the compressive sensing theory, Lustig et al., Sparse MRI: The application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine, 58, 1182-1195 (2007) (Lustig et al., 2007, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein) proposed effectively reconstructing MR images with only 20% sampling.
The improved results were obtained by having both a wavelet transform and a discrete gradient in the objective, which can be formulated as follows:
                              x          ^                =                  arg          ⁢                                          ⁢                                    min              x                        ⁢                          {                                                                                          1                      2                                        ⁢                                                                                                                    Rx                          -                          b                                                                                            2                                                        +                                      α                    ⁢                                                                                          x                                                                    TV                                                                      =                                  β                  ⁢                                                                                                          Φ                        x                                                                                    1                                                              }                                                          (        1        )            
where α and β are two positive parameters, b is the under-sampled measurements of k-space data, R is a partial Fourier transform and Φ is a wavelet transform. This reconstruction can be based on the fact that piecewise smooth MR images of organs can be sparsely represented by the wavelet basis and should have small total variations.
The total variation (“TV”) can be viewed discretely as:
                  x              TV    =            ∑      i        ⁢                  ∑        j            ⁢              (                                            (                                                ∇                  1                                ⁢                                  x                  ij                                            )                        2                    +                                    (                                                ∇                  2                                ⁢                                  x                  ij                                            )                        2                          )            
where ∇1 and ∇2 denote the forward finite difference operators on the first and second coordinates, respectively. Partial differential equation (“PDE”), (He, et al, MR image reconstruction by using the iterative refinement method and nonlinear inverse scale space methods. Technical Report UCLA CAM 06-35 (2006) (He et al., 2006, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein)) methods have been used to attack it. However, they are slow and impractical for real MR images. Computation need not become the bottleneck that can prevent such a model, as discussed in He et al., 2006, from being used in practical MR image reconstruction. Therefore, a problem to be addressed in compressed MR image reconstruction is to develop efficient algorithms to solve this problem with equivalent or better reconstruction accuracy.
Other methods tried to reconstruct compressed MR images by performing Lp-quasinorm (p=1) regularization optimization (See, e.g., Ye, et al., Projection reconstruction MR imaging using focuss, Magnetic Resonance in Medicine, 57, 764-775 (2007) (Ye et al., 2007, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein)); Chartrand, Exact reconstruction of sparse signals via nonconvex minimization, IEEE Signal Processing Letters, 14, 707-710 (2007) (Chartrand, 2007, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein); Chartrand, Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data, Proceedings of the Sixth IEEE International Conference on Symposium on Biomedical Imaging, IEEE Press, pp. 262-265 (2009) (Chartrand, 2009, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein). Although a bit higher compression ratio may be achieved, these non-convex methods do not always give global minima and are also relatively slow. Trzasko et al., Highly under-sampled magnetic resonance image reconstruction via homotropic l0-minimization, IEEE Transactions on Medical Imaging, 28, 106-121 (2009) (Trzasko et al., 2009, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein), discussing homotopic non-convex l0-minimization to reconstruct MR images. They created a gradual non-convex objective function which may allow global convergence with designed parameters. It was faster than those Lp-quasinorm regularization methods. However, it still needed 1-3 minutes to obtain reconstructions of 256×256 images in MAT-LAB on a 3 GHz desktop computer.
Recently, two fast methods were proposed to directly solve the above noted iterative refinement method and nonlinear inverse scale space MRI image reconstruction “Problem 1” solution model. In Ma, et al., An efficient algorithm for compressed MR imaging using total variation and wavelets, In IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2008, pp. 1-8 (2008) (Ma et al., 2008), an operator-splitting algorithm (TVCMRI) was proposed to solve the above noted MR image reconstruction problem. In Yang, et al. A fast alternating direction method for TV/1-/2 signal reconstruction from partial Fourier data, IEEE Journal of Selected Topics in Signal Processing, special Issue on compressive Sensing, 4, pp 288-297 (2010) (Yang et al., 2010, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein), a variable splitting method (RecPF) was proposed to solve the MR image reconstruction problem. Both can be used to replace iterative linear solvers with Fourier domain computations, which can gain substantial time savings. In MATLAB on a 3 GHz desktop computer, they can be used to obtain good reconstructions of 256×256 images in ten seconds or less. They are two of the fastest MR image reconstruction methods so far.
The iterative refinement method and nonlinear inverse scale space MRI image reconstruction solution model can be interpreted as a special case of general optimization problems consisting of a loss function and convex functions as priors. Two classes of algorithms to solve this generalized problem are operator-splitting algorithms and variable-splitting algorithms. Operator-splitting algorithms search for an x that makes the sum of the corresponding maximal-monotone operators equal to zero. These algorithms can use Forward-Backward schemes Gabay, Chapter IX applications of the method of multipliers to variational inequalities, Studies in Mathematics and its Applications 15, 299-331 (1983) (Gabay, 1983, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein); Combettes, et al., Signal recovery by proximal forward-backward splitting, SIAM Journal on Multiscale Modeling and Simulation 19, 1.107-1130 (2008), (Combettes, et al. 2008, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein); Tseng, A modified forward-backward splitting method for maximal monotone mappings, SIAM Journal on Control and Optimization, 38, 431-446, 2000 (Tseng, 2000, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein), Douglas-Rachford splitting schemes, Spingarn, Partial inverse of a monotone operator, Applied Mathematics and Optimization, 10, pp. 247-265 (1983), (Spingarn, 1983, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein) and projective splitting schemes Eckstein, et al., General projective splitting methods for sums of maximal monotone operators, SIAM Journal on Control Optimization, 48, 787-811 (Eckstein, et al., 2009, the disclosure of kyhich is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein).
The Iterative Shrinkage-Thresholding Algorithm (ISTA) and Fast ISTA (FISTA) Beck, et al., A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2, 183-202 (2009) (Beck, et al. 2009 II, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein) are two well known operator-splitting algorithms, which have been successfully used in signal processing. Beck, et al., Fast gradient-based algorithms for constrained total variation image demising and deblurring problems, IEEE Transaction on Image Processing 18, 2419-2434, 2009 (Beck, et al. 2009 I, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein), and Beck et al, 2009 II) and multi-task learning, Ji et al., An accelerated gradient method for trace norm minimization, in: ICML '09 Proceedings of the 26th Annual International Conference on Machine Learning, pp. 457-464 (2009), (Ji 2009).
Variable splitting algorithms, on the other hand, are based on combinations of alternating direction methods (ADM) under an augmented Lagrangian framework. It is firstly used to solve the PDE problem in Gabay, et al. A dual algorithm for the solution of nonlinear variational problems via finite-element approximations, Computers and Mathematics with Applications 2, 17-40 (1976) (Gabay, et al. 1976, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein); Glowinski, et al., Augmented Lagrangian and operator-splitting methods in nonlinear mechanics, Society for Industrial Mathematics, 9 (1989) (Glowinski, et al., 1989, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein). Tseng et al. extended it to solve variational inequality problems, Tseng, Applications of a splitting algorithm to decomposition in convex programming and variational inequalities, SIAM Journal on Control and Optimization, 29, 119-138, 1991, (Tseng, 1991, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein); He, et al. A new inexact alternating direction method for monotone variational inequalities. Mathematical Programming 92, 103-118 (He et al., 2002 II, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein). Wang et al., A new alternating minimization algorithm for total variation image reconstruction, SIAM Journal on Imaging Sciences, 1, pp. 248-272 (2008), (Wang et al. 2008, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein) showed that the ADMs are very efficient for solving TV regularization problems. They also outperform previous interior-point methods on some structured SDP problems, Malick, et al., Regularization methods for semidefinite programming, SIAM Journal on Optimization, 20, pp. 336-356, (2009), (Malick et al., 2009, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein).
Two variable-splitting algorithms, namely the Multiple Splitting Algorithm (MSA) and Fast MSA (FaMSA), have been recently proposed to efficiently solve the above noted iterative refinement method and nonlinear inverse scale space MRI image reconstruction solution model while all convex functions are assumed to be smooth Goldfarb, et al. Fast Multiple Splitting Algorithms for Convex Optimization, Technical Report, Department of IEOR, Columbia University, New York (2009) (Goldfarb, et al., 2009, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein).
However, all these above-mentioned algorithms cannot efficiently solve the above noted iterative refinement method and nonlinear inverse scale space MRI image reconstruction model problem with provable convergence complexity. Moreover, none of them can provide the complexity bounds of iterations for their problems, except ISTA/FISTA, as discussed in Beck, et al., 2009 II, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein, and MSA/FaMSA, as discussed in Goldfarb, et al. 2009, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein. Both ISTA and MSA are first order methods. Their complexity bounds are (1/ε) for ε-optimal solutions. Their fast versions, FISTA and FaMSA, have complexity bounds (1/√{square root over (ε)}), which are inspired by the seminal results of Nesterov and are optimal according to the conclusions of Nesterov, Nesterov, A method for solving the convex programming problem with convergence rate o(1/k2), Dokl. Akad. Nauk SSSR. 269, 543-547, 1983, (Nesterov, 1983, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein, Nesterov, Gradient methods for minimizing composite objective function, Technical report at www.ecore.beDPs/dp1191313936.pdf. Nesterov, 2007, the disclosure of which is hereby incorporated by reference in its entirety for all purposes, as if actually entirely reproduced herein).
However, both ISTA and FISTA are designed for simpler regularization problems and cannot be applied efficiently to the composite regularization problem iterative refinement method and nonlinear inverse scale space MRI image reconstruction model problem using both L1 and TV norm. While the MSA/FaMSA methods assume that all convex functions are smooth, it makes them unable to directly solve the iterative refinement method and nonlinear inverse scale space MRI image reconstruction model problem, as the non-smooth function must be smoothed first before applying them. Since the smooth parameters are related to ε, the FaMSA with complexity bound (1/√{square root over (ε)}) requires (1/ε) iterations to compute an ε-optimal solution, which means that it is not optimal for the iterative refinement method and nonlinear inverse scale space MRI image reconstruction model problem.
A formalism widely used in magnetic resonance imaging is k-space. Simply speaking, k-space is the temporary image space in which data from digitized MR signals are stored during data acquisition. When k-space is full (at the end of the scan), the data can be mathematically processed to produce a final image. Thus k-space holds raw data before reconstruction. The so-called k-space is in the spatial frequency domain. Thus if one defines kFE and kPE such that kFE= γGFEmΔt and kPE= γnΔGFEτ, where FE refers to frequency encoding. PE to phase encoding, Δt is the sampling time (the reciprocal of sampling frequency), τ is the duration of GPE, γ is the gyromagnetic ratio, in is the sample number in the FE direction and n is the sample number in the PE direction (also known as partition number), the 2D-Fourier Transform of this encoded signal results in a representation of the spin density distribution in two dimensions. Thus position (x, y) and spatial frequency (kFE, kFE) constitute a Fourier transform pair. The so-called k-space can have the same number of rows and columns as the final image. During the scan, k-space can be filled with raw data one line per TR (Repetition Time).
Although a strict mathematical proof does not exist and counterexamples can be provided, in most cases it is safe to say that data in the middle k-space (high frequency information) contains the signal to noise and contrast information for the image, while data around the outside of the image (low frequency information) contains all the information about the image resolution. This is the basis for advanced scanning techniques, such as the keyhole acquisition, in which a first complete k-space is acquired, and subsequent scans are performed by acquiring: just the central part of the k-space, In this way, different contrast images can be acquired without the need of running full scans.
A symmetry property exists in k-space, descending from the fact that the object imaged is a contrast-weighted proton density and thus a real quantity, relating the signal at two opposite locations in k-space:S(−kFE,−kPE)=S*(kFE,kPE)
where the star (*) denotes complex conjugation. Thus k-space information is somewhat redundant, and an image can be reconstructed using only one half of the k-space, either in the PE (Phase Encode) direction saving scan time (such a technique is known as half Fourier or half scan) or in the FE (Frequency Encode) direction, allowing for lower sampling frequencies and/or shorter echo times (such a technique is known as half echo). However, these techniques are approximate due to remaining phase errors in the MRI data which can rarely be completely controlled.